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ABSTRACT 

We present local simulations that verify the linear streaming instability that arises from aerodynamic 
coupling between solids and gas in protoplanctary disks. This robust instability creates enhancements 
in the particle density in order to tap the free energy of the relative drift between solids and gas, 
generated by the radial pressure gradient of the disk. We confirm the analytic growth rates found by 
Youdin & Goodman (2005) using grid hydrodynamics to simulate the gas and, alternatively, particle 
and grid representations of the solids. Since the analytic derivation approximates particles as a fluid, 
this work corroborates the streaming instability when solids are treated as particles. The idealized 
physical conditions - axisymmetry, uniform particle size, and the neglect of vertical stratification and 
collisions - provide a rigorous, well-defined test of any numerical algorithm for coupled particle-gas 
dynamics in protoplanctary disks. We describe a numerical particle-mesh implementation of the drag 
force, which is crucial for resolving the coupled oscillations. Finally we comment on the balance of 
energy and angular momentum in two-component disks with frictional coupling. A companion paper 
details the non-linear evolution of the streaming instability into saturated turbulence with dense 
particle clumps. 

Subject headings: diffusion — hydrodynamics — instabilities — planetary systems: protoplanetary 
disks — solar system: formation — turbulence 



1. INTRODUCTION 

Solid bodies in protoplanetary disks lose angular mo- 
mentum as they encounter the headwind of the pressure- 
supported gas disk. The subsequent radial drift is fastest 
for marginally coupled solids whose aerodynamic stop- 
ping times are comp arable to the local orbital time 
([Weidenschilling] I1977T ) . For standard disk models, cm- 
sized particles at 30 AU and m-sized bodies at 1 AU suffer 
drift times of only approximately 10 or 100 orbital peri- 
ods, respectively. Rapid infall imposes severe time-scale 
constraints on the growth into km-sized solid bodies, or 
planetesimals, by coagulation. Concerns about the in effi- 
ciency of sticking for macroscopic solids (|Benzll2000ft has 
also contributed to the concept of a "meter-size barrier" 
in planet formation (which should not be misinterpreted 
as implying that gro wth to meter sizes is easy, see e.g. 
iBlum fc Wu7rr][200l . 

Th e gravitational instabilit y hypothesis (jSafronovl 
Il969t iGoldreich fc Ward! I1973D postulates that a sedi- 
mcntcd mid-plane layer of small particles (perhaps mm- 
sized to match chondrules) will fragment directly into 
gravitationally bound planetesimals, avoiding the prob- 
lems with sticking efficiency and drift. However, disk tur- 
bulence acts to diffuse particles, inhibiting both their ver- 
tical settling to the midplane (jWeidenschilling fc Cuzzil 
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119931 : iDubrulle et~afl I1995Q and their abi lity to col- 
lapse into bound structures (jYoudinl l2005h . Even in 
a completely laminar disk, particle settling generates 
vertical shear in the orbital motion of the gas. This 
shear in turn triggers modified Kelvin-Helmholtz insta- 
bilities that de velop into turbulence, restricting further 
sedimentation dGoldreich fc Wardl[l973l : IWeidensc hilling 
119801 : ICuzzi et all 119931 ). This self-induced turbulence 
may not be able to prevent gravitational collapse if 
the soli ds-to-gas ratio is enhanced above Solar abun- 
dance s (lSekival[l998l:lYoudin fc Shul l200llGaraud fcLhl 
12004 IWeidenschilring||200(j| ). possibly due to photoevap- 
orati on of the gas-rich sur face layers of the stratified 
disk (jThroop fc Ballvl [2005) or to pile -ups of solids in 
the inner disk from particles that drift in m ore rapidly 
from the outer disk ([Youdin fc C hiang 2004) . Significant 
progress has been made in under standing the turbulence 
generated by particle settling (Ishitsu & Sekiya: 120031 : 
iGomez fc~ Ostrikcr 2005c [johansen et al.ll2006f ). However 
a simulation that incorporates the full 3D nature of these 
non-axisymmetric instabilities, with radial shear and the 
independent evolution of solids and gas, has not yet been 
performed. 

This paper addresses the related streaming instability 
([Youdin fc Goodmanll2005l hereafter referred to as YG) 
where vertical gravity is ignored in order to focus on a 
simpler manifestation of particle-gas coupling in Keple- 
rian disks. With no vertical shear present, the stream- 
ing instability is driven by the relative motion between 
solids and gas, which is predominantly radial for tightly 
coupled particles. The ultimate energy source, as with 
vertical shear instabilities, is the radial gas pressure gra- 



2 



Youdin & Johansen 



dient. Particle feedback on gas dynamics is important 
not just for establishing the (unstable) equilibrium, but 
also for generating escalating oscillations. Consequently, 
streaming instabilities trigger exponential growth of ar- 
bitrarily small particle density pert urbations, as shown 
by YG . The single- fluid treatment of lGoodman fc Pindorl 
( 2000) discovered a related boundary layer drag instabil- 
ity in stratified disks t hat c ould also concentrate par- 
ticles. Uohansen et al.l (|2006t ) found significant particle 
clumping in studies of Kelvin- Helmholz instabilities with 
particle feedback on the gas, which those authors hypoth- 
esized was a manifestation of non-linear streaming insta- 
bilities. The current study , including a companion paper 
(| Johansen fc Youdin]|2007| ). explores the consequences of 
streaming instabilities, and more generally the role of 
particle-gas coupling in protoplanetary disks. This paper 
demonstrates that our simulations faithfully reproduce 
the linear physics of the streaming instability, whether 
the solids are modeled as a fluid or Lagrangian particles. 

The paper is built up as follows. In we present 
the basic equations of our dynamical system and review 
the streaming instability. Section [3] describes the nu- 
merical methods, including the communication of drag 
forces between particles and a grid in i)3.2l Our main 
results, in 21 numerically confirm the linear stream- 
ing instability. In <J5] we analyze energy and angular 
momentum balance in a coupled two-fluid system. We 
discuss our results in The appendices contain an 
analysis of interpolation and assignment errors in differ- 
ent particle-mesh approaches to calculating drag forces 
(Appendix [X]) . a non-axisymmetric analytical problem 
used to test drag force assignment over shear-periodic 
boundaries ( Appendix IB]) , and a recipe to minimize Pois- 
sion noise in seeding linear particle density perturbations 
(Appendix [Cj . A companion paper, Johansen & Youdin 
(2007, hereafter referred to as JY), describes the full non- 
linear evolution of the streaming instability into turbu- 
lence. 



dw 3^ dw n ^ 

— — h (w ■ v)w llx— — = 2S2w v x 
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Transport terms on the left hand side of equations (fTJ) 
and @ include advection by the peculiar velocities, w, 
and by the Kepler shear, Vo. The right hand side of 
the equation of motion (Eq. [5]) contains Coriolis forces 
(as modified by Kepler shear) and drag acceleration rel- 
ative to the gas component with velocity u. We apply 
a linear drag force with constant friction time tj, valid 
for relat ively small particles in the Epstein or S tokes 
regimes (lAdachi et all 119761 IWeidenschillindH977h . Ep- 
stein's Law, Tj Ep ^ = p,R/(p g c s ) holds for particles of 
size R < A g , where A g w (r/AU) 2,75 cm is the mean free 
path of the gas molecules, c s is the gas sound speed, 
and p, is the internal density of rock/ice. Stokes' Law, 

Tf St ^ — T f i?/A g applies in the relatively narrow range 
A g ^2 5; Ag^K/cs, where vk = is the local Keplerian 
speed. Yet larger particles, R > \ g v\^/c s , trigger turbu- 
lent wakes with non-linear drag accelerations, which can 
not be modeled with a constant friction time. 3 Note 
that Stokes' Law is independent of gas density (since 
A g oc 1/pg). The dependence of Epstein's law on gas 
density fluctuations is neglected in our calculations as it 
is a small correction for low Mach number flow. 

The solid component does not feel a pressure gradient, 
neither from the gas, because the mass per solid parti- 
cle is so high, nor from interparticle collisions, because 
the number density is so low. Drag effects dominate colli- 
sional effects, since the collision time, t co \\ — p.i?/(p p c p ), 
is long with t co \\/Tf s» (p g / p p )(c s /c p ) 3> 1, even when the 
particle density is large, since the rms speed of particles, 
c p , is much smaller than the gas sound speed. 4 

For numerical work, we also use a Lagrangian descrip- 
tion of particle motion, see H3.ll 



2. STREAMING INSTABILITY: ANALYTICS 

2.1. Basic Equations 

We describe the local dynamics of the gas and 
solid component of a protopl anetary disk in the shear- 
ing sheet approximation fe.g. [Goldreich fc Lvnden-Belll 



1965| iGoldreich fc Tremaindll978ir The Cartesian coor- 
dinate frame corotates with the Keplerian frequency f2 
at an arbitrary orbital distance r from the central grav- 
ity source. The coordinate axes are oriented such that 
x points radially outwards, y points along the rotation 
direction of the disk, while z points vertically out of the 
disk, parallel to the Keplerian rotation vector ft. Our 
unstratified model omits vertical gravity. We measure 
all velocities relative to the linearized Keplerian shear 
flow in the rotating frame Vo = V y< oy = — (3/2)f2xy. 

2.1.1. Solids as a Fluid 

Analytic investigations are greatly simplified by treat- 
ing solid particles as a continuous fluid of density p p 
and velocity w, which evolve according to shearing sheet 
equations of continuity and motion 
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2.1.2. Gas Evolution 

The equations of continuity and motion for the gas 
read 
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Equation reduces to V • it = for an incompressible 
gas, as was considered in YG. The momentum equation 
contains advection and Coriolis forces as equation 
([2]). The main distinction between the two components 

3 The onset of turbulent wakes would be stalled to larger par- 
ticles if the relative velocity \u — w\ < In practice, how- 
ever, particles this large are weakly coupled and experience the 
full pressure-supported headwind. 

4 If the particles were large enough for non-linear turbulent drag, 
then collisional effects could only be safely neglected for p g > p p 
and/or if drift motions dominate particle random motions. With 

fturb) ^ p m R/( p \ u - w \) t then t co il/T- f (tUrb) ~ (Pg/Pp)([«— H/°p)- 
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Fig. 1. — Linear growth rate s of the streaming instability vs. radial and vertical wavenumbers for a friction time of r s = 1.0 (upper row) 
and T s = 0.1 (lower row). Three values of the solids-to-gas density ratio, e = 0.2, 1.0, 3.0, are considered along the columns. Contours label 
log 10 (s/i7), darker shading corresponds to faster growth rates, while the dotted regions contain only damped modes. 



is that gas is effected by pressure gradients. We include 
both local pressure gradients from isothermal gas den- 
sity fluctuations and a constant acceleration by a global 
radial pressure gradient, 5 dP/dr, expressed using the di- 
mensionless measure of sub-Keplerian rotation 



n 



dP/dr 
2p g Q 2 r 



(5) 



The feedback of the linear drag force scales with the den- 
sity ratio of particles to gas, 



e = Pp/Ps > 

which ensures that total momentum is conserved. 



(6) 



2.2. 



Equilibrium State 

Equilibrium solutions to the mutually coupled equa- 
tions ^ and (0| were obtained bv lNakagawa et aLl (1986, 
hereafter referred to as NSH) for local and linear dynam- 
ics. The in-plane deviations from Keplerian rotation are 



2eT s 
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rV v K ■ 



(7) 



5 While seemingly paradoxical, it is consistent with the local 
approximation to include global pressure gradients to linear order 
while treating the background gas density as constant. 
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(10) 



The dimensionless stopping time, r s = f?Tf, is a conve- 
nient measure of coupling strength, since marginal cou- 
pling, r s = 1, famously maximizes the radial drift speed 
of an isolated particle. Velocities scale with the sub- 
Keplerian velocity, t]Vk, where = fir. 6 The azimuthal 
velocities are factored into the center-of-mass motion, 



^(com) _ Pg U V + Pp W V 

V Pp + Ps 



l + e 



(11) 



and order r s 2 drift motions (see YG for details). 
Vertical gradients in the solids-to-gas ratio e give gra- 

(for t s <1) that trigger the 



(com) 



dients in V % 

settling-induced Kelvin-Helmholz instabilities discussed 
in the introduction. As in YG, we also neglect vertical 

6 The radial dependence of equilibrium quantities, such as rjvj^ 
and T s , is ignored in the local approximation since the effect of the 
Keplerian shear profile is already included (to linear order) in the 
equations of motion. 
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gravity in the present work in order to allow for a lam- 
inar equilibrium state. With vertical gravity, any ini- 
tial condition must be time-dependent (due to vertical 
settling) and/or turbulent (to halt the settling). Fur- 
thermore, in stratified disks, drift speeds (and even di- 
rections) vary with height above the midplane, since r s 
rises with decreasing gas density and since the radial gas 
ressure gradient can reverse away from the mid-plane 
Takeuchi fc Lin! [20021 ). This is particularly relevant for 
small grains that remain above the midplane for many 
orbital times. The severity of the unstratified approxi- 
mation is justified by the insights gained from an initially 
simple, well-defined problem that rapidly turns complex. 



2.3. Streaming Instability 

The streaming motion of solid particles through gas 
presents a source of free energy that is driven by pres- 
sure gradients and mediated by drag and Coriolis forces. 
YG showed, by linearly perturbing equations |(5J) and ^ 
about the equilibrium state given by equations ((7J)— (JTUJ) , 
that this streaming robustly triggers instability in proto- 
planetary disks. The instability provides a novel mecha- 
nism to generate growing particle density perturbations 
in a moderately dense mid-plane layer of macroscopic 
particles, while smaller particles (r s <C 1) with poor drag 
feedback (e <C 1) will give rise to only very low, sub- 
dynamical growth rates. 

The YG analysis and the linear test simulations in this 
paper are "2.5-D", i.e. all three components of velocity 
fluctuations are considered, 7 but perturbations are ax- 
isymmetric and characterized by the radial and vertical 
wavenumbers, k x and k z . The growth rates for several 
choices of r s and e (which henceforth indicates the aver- 
age value of Pp/pg in the background state, unless oth- 
erwise noted) are shown in Fig. Q] as a function of the 
dimensionless wavenumbers K x = k x r\r and K z = k z r\r. 

Since particles only affect gas dynamics via drag feed- 
back, growth rates increase for larger e, while the relevant 
length scales shrink, most likely because the response 
time-scale of the gas speeds up as Tf/e. Fig. [2] shows 
these trends, along with the particularly sharp increase of 
s across e = 1 for tightly coupled particles with r s = 0.1. 
The crucial physical distinction for marginal coupling 
(for which the same sharp increase is not present) may 
be that for r B « 1, azimuthal drift (of order t s 2 ) is no 
longer negligible compared to radial drift (of order t s ). 
For a more technical difference, note the gray curves in 
Fig. [21 which show that the phase speed of waves changes 
sign near e as 1. YG noted that the phase speed tends to 
track the component with the fastest radial drift - solids 
for e < 1 and gas for e > 1. Curiously at r s = 1 the 
transition is delayed to e ~ 2. As r s decreases the switch 
in phase speeds gets closer to e = 1, coinciding with the 
rise in growth rates across e = 1 becoming steeper and of 
larger amplitude (see also Fig. 3 of YG for the r s = 0.01 
case). 

The trend with r s is complicated as well. In the gas- 
dominated regime (e < 1) growth rates show the ex- 
pected rise toward the r s « 1 "sweetspot" : streaming 
motions are large yet particles still respond effectively to 
the gas. The situation reverses when particles dominate 

7 And all three components are necessary for axisymmetric in- 
stability (YG). 
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Fig. 2. — Peak growth rate, s, of the streaming instability and 
fastest growing radial wavenumber, k x , versus the solids-to-gas 
density ratio e = p p /p g for a friction time of r s = 1.0 (solid line) 
and t b =0.1 (dashed line). Growth becomes faster and occurs at 
smaller scales for increasing e, with a particularly sharp increase 
in s across e = 1 for tightly coupled particles with t s = 0.1. Gray 
curves in lower plot (associated with gray axis on right) show the 
radial phase speed of waves. The sharp dips near e si 1—2 indicate 
a sign change for the wave speed: inward when gas dominates and 
outward when particles dominate. 



(e > 1), with growth rates that are actually faster for 
tighter coupling, but at smaller length scales. 

Returning for a moment to Fig. [H it is also evident 
that growth does not peak at a single pair of wavenum- 
bers. The fastest growing K x can be determined, with 
only damped modes for sufficiently large K Xl but growth 
remains flat for large K z (indeed the curves of Fig. [2] are 
calculated in the limit K z /K x ^ 1). A physical expla- 
nation for the difference between large or small K z /K x 
follows. The (near) incompressibility of the gas imposes 
a ratio Imz/m^I — \K X /K Z \. With K z 3> K x , velocity 
vectors arc nearly parallel to the x — y plane with neg- 
ligible vertical velocities (just enough to maintain gas 
incompressibility). Since the balance of Coriolis forces is 
maintained in thin vertical sheets, instability persists to 
large K z . On the other hand, large K x /K z shrinks u x 
and destroys the necessary balance of Coriolis forces. 

The linear growth regime is surprisingly complex, con- 
sidering the simplicity of the physical system. Toy mod- 
els to explain the mechanism have unfortunately fallen 
short of capturing the essence of the instability. For in- 
stance, one might suspect that, since streaming insta- 
bilities involve particle density enhancements, they arise 
because radial drift slows in overdense regions [see equa- 
tion ©] leading to local traffic jams. 8 This effect, while 
relevant, does not explain linear growth of infinitesimal 
perturbations. To see this, consider the axisymmetric 
evolution of particle density that follows from the equi- 
librium drift speed [equation l[9])] and continuity [equa- 
tion (JTJ], which we express for simplicity in terms of a 

8 Not to be confused with t he global pile-ups t hat arise 
for Epstein drag either without ((Ypudin & Shu 2002) or with 
fYoud in & Chiang 2004) drag feedback. 
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TABLE 1 

Test Mode Eigensystems 































u x 


Uy 


u z 


Pg 




Wy 




UJ 


linA: t b 


= 0.1, £ 


= 3.0 


-0.1691398 


+0.1336704 


+0.1691389 


+0.0000224 


-0.1398623 


+0.1305628 


+0.1639549 


-0.3480127 


(K x = 


= 30, K z 


= 30) 


+0.0361553i 


+0.0591695i 


-0.0361555i 


+0.0000212i 


+0.0372951i 


+0.0640574i 


-0.0233277i 


+0.4190204i 


linB: r s 


= 0.1, e 


= 0.2 


-0.0174121 


+0.2767976 


+0.0174130 


-0.0000067 


+0.0462916 


+0.2739304 


+0.0083263 


+0.4998786 


(K x = 


= 6, K z -- 


= 6) 


-0.2770347i 


-0.0187568i 


+0.2770423i 


-0.0000691i 


-0.2743072i 


+0.0039293i 


+0.2768866i 


+0.0154764i 



Note. — Frequency ui is normalized to Q, velocities are normalized to and densities to the average value for particles or gas respectively. 
All eigenvalue coefficients are relative to the particle density perturbation, which should be set to p p < 1 for the evolution of the mode to be 
linear. We used p p = 10 -6 to normalize the eigenvector. The (tiny) effect of compressibility is included in the coefficients with r]v^/c B = 0.05. 
The growth rate s is the imaginary part of uj. 



variable (only for now) e = p p (x,t)/ p g .o = eo + e'(x,t) as 



dt 



dx 



(l + e) 2 +' 



(12) 



Linearizing about e' -C £o clearly gives stable wave prop- 
agation at the drift speed. Non-linear perturbations in 
equation (|12j) will steepen a particle density wave, with 
no amplitude grow th (readily shown by the method of 
characteristics, see Shu 1992). Even if the traffic jam 
concept fails to explain the linear growth of the stream- 
ing instability, it may be used to explain the non-linear 
clumping seen in JY (see also ^5.11 in this paper). 

We find in JY that non-linear states also show remark- 
able diversity with friction time and solids-to-gas ratio. 
We must, however, first ensure that the numerical algo- 
rithms can capture and confirm the linear growth phase. 

2.3.1. Eigenvectors and Vertical Standing Waves 

To test the growth rates of Fig. [T] computationally, 
the eigenvectors, i.e. relative amplitudes and phases of 
the density and velocity perturbations, must be carefully 
seeded for a specific choice of parameters r s ,e, K x , K z . 
The perturbation in each dynamical variable / can be 
written in terms of its complex amplitude / (a compo- 
nent of the full eigenvector) as fix, z) = $l{f exp{i(k x x + 
k z z — cot)]}, where to = lus% + is is the complex eigen- 
value containing the wave frequency wjf and the growth 
rate s. We choose to eliminate the superfluous vertical 
phase speed by superposing pairs of modes with vertical 
wavenumbers k z and — k z , respectively. Under a vertical 
parity transformation the vertical velocity amplitudes are 
odd, while all others are even. The superposition yields 

f e (x,z) = [3?(/) cos(k x x - ujftt) - 

3(/) sin(fc x x — LJsRt)] cos(k z z) exp(st) , (13) 
fo(x, z) = -[SR(/) sin(k x x - uj^t) + 

9f(/) cos(k x x — ujfRt)] sin(k z z) exp(si) , (14) 

for even (e) and odd (o) dynamical variables, respec- 
tively, which are now clearly standing waves in z. 

Table [1] lists eigenvalues and eigenvectors for the cases 
we will test numerically in 21 The calculation is similar 
to that of YG except gas compressibility was added so 
that a gas density perturbation can be included in the 
numerical calculations. The effect of the gas compress- 
ibility is otherwise negligible for t/vk/cs ~ c s /vk *C 1 
(the reason it was neglected in YG), affecting eigenval- 
ues and eigenvectors in the 5th digit for our choice of 
t]Vk/c s — 0.05. We also checked that the sound waves in- 
troduced by gas compression are rapidly damped. Note 



that Table Q] shows the gas density (and thus pressure) 
perturbations are out of phase (by ~ 90° and ~ 180° 
for A and B, respectively) with the particle density per- 
turbation. Thus solids are not merely collecting in pres- 
sure maxima, as occurs in gas density structures that are 
steady in time. 

3. NUMERICAL METHODS 

As a numerical solver we use the Pencil Code. 9 This 
is a modular finite difference code that uses 6th order 
symmetric spatial derivativ es and a 3rd order Runge- 
Kutta time integration (see lBrandenburgll2003L for de- 
tails). A module already exists for solving the equa- 
tion of motion of a dust fluid th at interacts with the 
main gas fluid through drag force (jJohansen et alJfeOO^ 
lJohansen"fc Klahr 2005). The basic dynamical equations 
in the (here unstratified) shearing sheet are equations §3§ 
and ([4} for the gas and equations (Qj) and ([2]) for the 
solids. This equation set is stabilized by adding small 
diffusive terms to the equation of motion and by upwind- 
ing the ad vection term in the cont i nuity equations (fo r 
details, see lJohansen fc Klahrll2005 iDobler et al.ll2006h . 
Treating particles as a fluid facilitates analytic calcula- 
tions and is significantly cheaper for numerical simula- 
tions, but is not always the desired approach. 

3.1. Solids as Particles 

Using Lagrangian particles provides a more realistic 
description of the dynamics of the solids, and there are 
two main reasons to justify the additional effort. 10 First, 
particles at a given position need not have a single well- 
defined velocity as the fluid approximation assumes, i.e. 
trajectories can cross. This concern is particularly valid 
for marginal and looser coupling. Second, and more se- 
riously, the fluid treatment cannot capture large density 
gradients, especially since the "sound speed" of the pres- 
sureless fluid is zero. Stabilization of steep density gra- 
dients would require a large artificial viscosity that com- 
promises the dynamics. Thus a Lagrangian treatment 
of the solids is necessary for the non-linear simulations 
of JY which generate large particles overdensities. Since 
the analysis of YG describes solids as a fluid, we must 
demonstrate that the instability does not depend cru- 
cially on this assumption. 

When treating solids as numerical particles, or rather 
as superparticles since each numerical particle effectively 

9 The code is publicly available at 

http : //www. nordita. dk/data/brandenb/pencil- code/ 

iu See Garaud et al. (2004) for a thorough analysis of the validity 
of fluid descriptions of particle motion subject to gas drag 
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represents a huge number of individual solids, each par- 
ticle i has a position x^ and a velocity relative to 
the Keplerian s hear. Particle motions are governed by 
Hill's equations (j Wisdom fc Tremaindll988l ) 



dw« 
dt 



dccM 
dt 



1 



=t> « 



2 



/*>-u(o;W) , (15) 



(16) 



here including drag force and expressed in a form to ap- 
pear as the Lagrangian equivalent to equation For 
axisymmetric simulations in the radial-vertical plane, the 
evolution of v y (t) is included but the azimuthal compo- 



nent of equation (fl6|) is irrelevant, effectively replaced 
by dy^'/dt = since that dimension that is not present. 
The interpolation of gas velocities at the particle posi- 
tions, u(ccW), is addressed in the next section. 

3.2. Drag Force Calculation 

The computation of drag forces between Lagrangian 
particles and an Eulerian grid requires some care to avoid 
spurious accelerations and to ensure momentum conser- 
vation. Small errors in the gas velocity are dangerously 
amplified by the subtraction of highly correlated particle 
velocities. Our drag force algorithm involves three steps: 

1. Interpolating gas velocities at particle positions 

2. Calculating the drag force on particles 

3. Assigning the back-reaction force to the gas from 
particles in nearby cells 

For the first step, interpolation, we begin with gas ve- 
locities, , defined on a uniform grid where the index j 
labels the cells centered on positions x^ . We interpolate 
to the particle positions, x^ l \ using a weight function, 
Wi, as 

it(a;M ) =V Wi(x® -x^)u {j) 



3 



(17) 



The weight function is normalized as . Wi(x^ — 
xU>) = 1 for any x^\ and has non-zero contributions 
only from the cells in the immediate vicinity of x^ 3 \ 

The second step, calculating the drag acceleration on 
particle i, 

Tf 

is trivial once the relevant quantities are defined, but this 
is the step that amplifies interpolation errors in u(xW), 
because of strong coupling to particle velocities, a prob- 
lem that worsens for smaller Tf. We note that other 
choices of the drag law (e.g. non-linear in the velocity 
or including gas density fluctuations in Epstein drag) 
would be simple to implement by interpolating the rele- 
vant grid-based quantities as in equation (|17p . 

Finally, we calculate the back-reaction drag force, / g , 
on the gas in cell j. Assigning particle velocities to a 
mesh risks violating momentum conservation. Instead 



we follow the suggestion of Jim Stone (personal commu- 
nication) and use Newton's third law to directly assign 
the force on the particles back to the gas, 



ptf) 



Ps V cett i 



W A (x 



(i) 



)f ( v> , (19) 



where m p is the mass of a particle (if not uniform it 
would be inside the sum), and V ce \i is the volume of a 
grid cell. The assignment function Wa obeys the same 
conditions as Wi, so that only particles in a given cell or 
its nearby neighbors contribute to the sum. Global mo- 
mentum conservation follows trivially from summation 
of equation p^|) . 



(20) 



with no reference to the drag law, the interpolation func- 
tion, or any properties of the assignment function except 
normalization. Thus unlike particle-mesh calculations 
with interacting particles (e.g. by self-gravity), we are 
flexible to choose W\ and Wa independently, without vi- 
olating momentum conservation. Nevertheless, choosing 
Wa — W\ is safest since drag forces from gas to particles 
- and vice-versa - are smoothed symmetrically. 

We opted for second order interpolation and assign- 
ment methods, either quadratic spline or quadratic poly- 
nomial, which use three grid cells in each dimension, for 
a total of 9 (27) for 2-D (3-D) simulations, respectively. 
This gave considerable improvement over lower order bi- 
linear interpolation (but at a computational cost - the 
drag force calculations dominate the wall time in our sim- 
ulations with high order interpolation and assignment). 
The details and errors associated with the interpolation 
schemes are described in Appendix |Al The quadratic 
spline assignment /interpolation method is often referred 
to as the Triangular Shap ed Cloud scheme (TSC, see 
iHocknev fc Eastwoodlli981l ). 

3.2.1. Boundary Conditions for the Drag Force 

Our implementation of periodic boun dary conditions, 
and use of higher (than zeroth, as in iJohansen et al.l 
120061 ) order assignment schemes, causes particles near 
grid edges to exert drag forces on mesh points across the 
boundaries. In non-axisymmetric simulations (such as 
the 3-D simulations that we present in JY) the radial di- 
rection is shear-periodic so that two connected points 
at the inner and outer radial boundary are Ay(t) = 
mod[(3/2)f2L x t, L y ] apart in the azimuthal direction. 
Techniques for implementing radial boundary condition s 
in the shearing box are well-known (|Hawlev et al.|[l995j) . 
Fluid variables in zones on one radial boundary are 
copied to ghost zones adjacent to the opposite boundary 
and shifted azimuthally. Then differences across bound- 
aries are performed, i.e. "copy, shift, and difference." 

The implementation of shear periodic boundary con- 
ditions for drag forces on the gas is a subtly different 
"assign, shift, and add" procedure, as sketched in Fig. [3l 
First we assign the (appropriate fraction of) drag accel- 
erations from particles in boundary zones to gas in the 
ghost zones. Then we shift the accelerations on the ra- 
dial ghost zones in the y-direction, the inner by —Ay(t), 
and the outer by +Ay(t). Finally these shifted accel- 
erations are added (or folded) to the first real zone on 
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Fig. 3. — A sketch of the shear-periodic radial (x) boundary con- 
dition for the assignment of drag forces from a particle to the gas. 
The dot represents a particle near the boundary and crosses indi- 
cate the (centers of) gas cells that receive a drag acceleration with 
the second order TSC assignment scheme (grayscale of crosses in- 
dicates rough weight of drag force received by gas in each cell) . We 
illustrate an example with 4 processors in the y-direction (labeled 
ipy). The periodic direction is indicated by solid diagonal lines. 
The drag force assigned to ghost cells across the boundary (circled 
on left) is shifted in Fourier space and then added as an acceleration 
on the physical grid cells at the outer boundary. Note that in prac- 
tice (a) the drag force from an individual particle influences more 
than three grid cells across the boundary, since displacements are 
not integer multiples of the grid spacing and (b) drag forces from 
all particles on a ghost zone are added before Fourier shifting. 

the opposite side of the mesh. We interpolate (since the 
ghost zones do not slide by integer numbers of grid cells) 
by applying the azimuthal shift in Fourier space. Fourier 
interpolation has the advantage over high order polyno- 
mial interpolation that the function and all its derivatives 
are continuous. A numerical test of the radial boundary 
condition with shearing waves is described in Appendix 

m 

4. NUMERICAL TESTS OF LINEAR GROWTH 

We now present measurements of linear growth rates 
of the streaming instability from numerical simulations. 
These results confirm the capabilities of our code and 
verify the authenticity of this fundamental instability 
not yet explicitly established for a particle-based treat- 
ment of solids. Our efforts in reproducing growth rates 
to a satisfactory accuracy were useful in developing our 
numerical implementation of drag forces. We hope that 
others who simulate coupled particle-gas disks will con- 
duct similar dynamical tests of the simplest (identified) 
aerodynamic drag instability. 11 

We choose two different test problems: an eigenvector 
for r s = 0.1, e = 3.0, K x = K z = 30 (run linA), which 
grows rapidly with s/Q = 0.41902, and an eigenvector 
for r s = 0.1, e = 0.2, K x = K z = 6 (run linB) that 
grows more slowly with s/Q — 0.01548 and hence is more 
numerically demanding. The total initial velocities are 

11 AY will provide initial conditions (eigenvectors) by e-mail 
request. 



the sum of the equilibrium drift solutions of equations 
(|7])- (fT0|) . and the vertically standing wave of equations 
(fT5|) - ([Li)) with eigenvectors from Table Q] The initial 
amplitude of the particle density was set to 10 -6 in all 
cases to ensure linearity. 

4.1. Growth for Solids as a Fluid 

The measured growth rate when particles are treated 
as a fluid is shown with a solid black line in Figs. [4] and 
[5] (the top and bottom plots are identical for the two- 
fluid case). The eight panels show the growth rate of 
the velocity and density of the gas (top row) and of the 
solids (bottom row) as a function of the number of grid 
points per wavelength. We have varied the resolution be- 
tween 3 and 64 grid points per wavelength for the fluid 
treatment of solids and between 8 and 64 grid points per 
wavelength for the particle treatment. The growth rates 
are obtained by spatially Fourier transforming the 8 dy- 
namical variables at 10 fixed times over At = 0.2i7 — 1 and 
measuring the amplitude growth of the relevant Fourier 
mode. There is generally an excellent agreement between 
the measured growth rates when the solids are treated as 
a fluid and the analytical values down to 4 grid points 
per wavelength, except for the gas density which shows 
some variation from the analytical value for crude reso- 
lutions. This disagreement is not surprising since small 
errors in the cancellation of du x /dx and du z /dz for the 
nearly incompressible gas give spurious growth to the gas 
density according to the linearized continuity equation 
d\np' g /dt = — V • u 1 . While the gas density perturba- 
tions are too small to affect the drag force, they also 
cause the pressure perturbations which are significant. 
Fortuitously, the errors in the gas density (for crude res- 
olutions) do not affect the other dynamical variables. It 
may help that spurious sound waves damp rapidly (in a 
stopping time). 

4.2. Growth for Solids as Particles 

Reproducing analytic growth rates using a particle rep- 
resentation of the solids is significantly more difficult 
than in the two-fluid case. Poisson fluctuations from 
undersampling and truncation errors in the drag force 
calculation cause numerical discrepancies. Section 13.21 
and Appendix [A] describe the algorithms for computing 
drag forces and the errors associated with interpolation 
and assignment. 

4.2.1. Cold Start Initialization 

To avoid shot noise in seeding linear particle density 
perturbations we use a "cold start" algorithm (described 
in detail in Appendix |C|) for the initial particle positions. 
First we place all particles on a uniform grid. Then we 
apply a small, spatially periodic shift to their positions. 
This seeds the desired mode with minimal noise leaked 
to other wavelengths. We experimented with different 
numbers of particles: 25 particles per grid cell to match 
the non-linear runs of JY, and 1 particle per grid cell as 
a test. 

With the cold start to eliminate noise and the TSC as- 
signment scheme to smoothly distribute a particle's in- 
fluence over the nearest three grid cells per dimension, 
communicating initial density perturbations of infinites- 
imal amplitude with only a few particles is trivial. Fig. [6] 
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Fig. 4. — Measured growth rate of a seeded mode with K x = K z = 30 and r s = 0.1, e = 3.0 as a function of the number of grid-points 
per wavelength, shown for 1 and 25 particles per grid cell in the top and bottom plots, respectively. The behavior of each dynamical 
variable is shown separately. The analytical growth rate, s = 0.419031?, is indicated with a gray line. The fluid treatment (solid black 
line) gives excellent agreement with the analytical growth rate down to 4 grid points per wavelength, whereas 16 grid points is needed for 
the regular TSC scheme (dash-dotted line). Applying Fourier sharpening to the initial condition gives some improvement (dashed line). 
Replacing spline interpolation with polynomial interpolation (dotted line) gives better growth rates, but polynomial interpolation has the 
disadvantage of being discontinuous over cell interfaces. Increasing the number of particles per grid cell from 1 to 25 has minimal influence 
on the linear growth. 

demonstrates the algorithm effectiveness with the near 
perfect replication of a 1-D particle density perturbation 
of amplitude 10~ 6 with only 32 grid cells and one par- 
ticle per cell. This is nothing more (or less) than the 
miracle of continuous numbers. The use of many parti- 
cles per grid cell is still necessary to get good statistics 
in non-linear simulations. 



4.2.2. Results 

The growth rates with solids as particles are shown 
(together with the two-fluid results) in Figs. [4] and [5] as 
a function of spatial resolution. The top and bottom 
plots in each figure are for 1 and 25 particles per grid 
cell, respectively. Particle number makes little difference 
for the agreement with linear theory, although additional 
particles give some improvement, notably for the growth 
rate of p p in Fig. O 

While all runs use the TSC scheme to assign drag forces 



to the gas, three different techniques were tested for the 
interpolation of gas velocities to particle positions: (1) 
quadratic spline interpolation, (2) quadratic spline inter- 
polation with an initial Fourier sharpening of the gas ve- 
locity field, and (3) quadratic polynomial interpolation. 
Errors in gas velocity interpolation are the most danger- 
ous since they are amplified in the force calculation by 
subtracting a particle velocity that is highly correlated 
with the gas flow. 

The first technique, quadratic spline interpolation, uses 
the same weight function as TSC assignment and gives 
smooth interpolates with a reduced fluctuation ampli- 
tude. The dash-dotted lines in Figs. H] and [5] show that 
this technique accurately reproduces the growth of p p . 
The results for the other variables are poor for resolu- 
tions of less than 16 grid points per wave length. This is 
a result of spurious drag forces generated because inter- 
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Fig. 5. — Similar to Fig. |4]but for a mode with K x = K z = 6 and r s = 0.1, e = 0.2 that has an analytical growth rate of s = 0.01548J7. 
The agreement with the analytical growth shows a comparable resolution dependance to Fig. [4] but here the increase to 25 particles per 
grid cell shows better agreement for p p . 



polation reduces gas fluctuation amplitudes. 12 

The second interpolation technique (shown with 
dashed lines in Figs. Q]and still uses quadratic splines, 
but sharpens the initial gas velocities to correct the drag 
force. The amplitude of the Fourier modes u are in- 
creased by the precise amount, [1 — A 2 (fc 2 + fc 2 )/8] _1 , 
that interpolation reduces them (see Appendix IX| . The 
sharpened TSC scheme gives much better growth rates, 
but still not as good as the two-fluid results. In a non- 
linear simulation with an evolving power spectrum, one 
could sharpen u with a pair of Fourier transforms at 
each time-step, but this was deemed too computation- 
ally costly. By getting improved results with only the 
initial condition sharpened, we show that growth rate 
discrepancies with spline interpolation are largely due to 
differences between numerical (discretized) and analytic 
eigenvectors that should not compromise the non-linear 
simulations. 

12 This is why, for unsharpened spline interpolation, growth rates 
are too large for u (gas is accelerated toward the unsmoothed am- 
plitude by particles) and too small for w (particles are decelerated 
by the lowered gas amplitudes). 



The third approach (shown with dotted lines in Figs. [4] 
and [5|) opts for precise quadratic polynomial interpola- 
tion instead of smoother splines. The resulting growth 
rates are comparable, or slightly better than, the sharp- 
ened splines. Despite the simplicity and good results 
obtained with this technique, we did not use it in the 
non-linear runs. Discontinuities in the interpolates at 
cell boundaries would add noise by leaking power to the 
grid scale. Since the errors of TSC are well-behaved 
(spatially smooth across a grid cell, declining with in- 
creasing resolution, and leaving particle density growth 
unaffected even at low resolution), we used spline in- 
terpolation in the non-linear runs. We also prefer the 
symmetry of using the same weight functions for inter- 
polation (quadratic spline) and assignment (TSC). 

Overall, numerical growth rates with solids treated 
as particles agree well with linear theory down to 16 
grid points per wavelength, although the particle den- 
sity grows at the correct rate even at 8 grid points per 
wavelength. Anomalies, particularly in the gas density, 
suggest that sound waves are being triggered due to in- 
terpolation errors, but these spurious motions damp and 
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Fig. 6. — A sinusoidal particle density perturbation of amplitude 
10 -6 as generated by the shifting algorithm of Appendix ICl with 
32 particles — only one per grid cell! The crosses (connected by the 
solid line) plot the TSC assignment of particle (over)density to the 
grid cells. Dots indicate the positions of the particles, but the shift 
is imperceptibly small. 

do not impede the expected growth of particle density 
perturbations. 

5. ENERGY AND ANGULAR MOMENTUM BALANCE 

This section provides brief overviews of energy and an- 
gular momentum in coupled particle-gas disks in order 
to provide a point of reference to more familiar dynam- 
ical systems, and because it will help us interpret the 
non-linear results of JY. We denote C = p g u y + p p w y as 
the total angular momentum density of solids and gas, 
ignoring the radius factor that is constant in the local 
approximation. The azimuthal components of equations 
© and © give 



d_ 

~dt 



— Ox-— 
2 dy 



C + V • T L 



-To 



dP 

dy 



(21) 



The terms on the left hand side relate local changes in C 
to the transport of C by the Keplerian flow and to the 
angular momentum flux T c = p g u y u + p p w y w. We do 
not call this flux a Reynolds stress because the velocities 
u and w have not been decomposed into fluctuations 
about their mean. The NSH equilibrium of equations 
(|7|)- (|10p transports angular momentum radially inwards, 



,u x u y + p p w x w y 



-2r s 3 p p 



rjv K 



(l + 6) 2 +T s 2 



(22) 

(23) 



a consequence of the slower rotation of the outgoing gas 
relative to the faster rotation of the incoming particles. 
This differs from the usual outw ard transport of angu- 
lar m omentum in accretion disks (|Lvn dcn-Bcll fc Pringld 
11974T ) , because the driving agent is not orbital shear, but 
the radial pressure gradient. 

The terms on the right hand side of equation (|21[) rep- 
resent sources or sinks of angular momentum: the ra- 
dial mass flux, T P}X = p g u x + p p w x , and azimuthal pres- 
sure gradients, where P is promoted to denote the to- 
tal gas pressure (background and perturbations) in this 



section. Equation (|2ip proves that axisymmetric equi- 
librium solutions cannot transport mass radially in the 
local model, a condition obeyed by equations ([7]) and ©. 
Note that equation (|2Tj) does not explicitly include drag 
forces, which transfer momentum between gas and solids, 
but (of course) do not dissipate C. 

The evolution of kinetic energy density £ = (p g |n| 2 + 
p p \w\ 2 )/2 is found by summing the dot products of p g u 
with equation (|4|) and p p w with equation (|2J) to give 



d_ 

dt 



3 o d 
2 dy 



£ + V-F £ =£ dl&g -u-VP+-nF c , x , 

( 24 ) 

where the energy flux, Ts = p g \u\ u + p p \w\ w, trans- 
ports energy radially inward (outward) when gas (parti- 
cles) dominate the mass, respectively. 13 The sources and 
sinks on the right hand side include the energy lost to 
drag dissipation, 



£a 



rag - 



-P P \w 
4(1 + 



[(1 



u\ 2 / Ti (25) 

\2 3 

" (r,v K ) 2 p p f2, (26) 



r 2 ) 2 



where the second equality applies to the NSH equilib- 
rium. A simple estimate of the effective temperature 
produced when the dissipated kinetic energy is released 
as thermal heat gives 

Tdrag < [£ p Q{riv K ) 2 la S v\ 1/4 ~ 30(r/AU)- 3 / 4 K (27) 

as an upper limit for the case of marginal coupling 
and e <^ 1, where S p ~ P P H P is the surface density 
of the solid component and H p is the scale height of 
the sublayer of solids. The above temperature limit is 
significantly colder than ev en passively irradiated disks 
([Chiang fc Goldreich1ll997h . a comforting fact for SED 
modelers. 

The £ W ork = —u ■ VP term represents energy gained 
from the work done by the total pressure forces. The 
equilibrium value of 



fwork = ~u x (dP/dr) = 



4r s 



(1 



-(t^k) 2 P p ^ (28) 



shows that |£ WO rk| > |£drag|, i-e. more energy is put into 
the system by pressure work than removed by drag. The 
final term, 



£ c ee (3/2)0^ = - 



3Tf 



[(l + c) 2 + T f 



2(vvk) 2 P p O , 



( 29 ) 

is well known in studies of viscous or collisional disks 
as the heat gen erated by the outward transport of angu- 
lar m omentum (jShu fc Stewartj fl985: Lit hwick fc Chiangl 
I2006h . However, in our case angular momentum trans- 
port is reversed according to equation (|2"3"]) and provides a 
sink of kinetic energy. The phenomenon of "backwards" 
angular momentum transport, and the dynamical cooling 
it provides, has been famously offered as an explanation 
for th e sharp edges of planetary rings (Bor deries et al.l 
[1981 . 

Equations ([26j) , p8]). and ([29]) verify that the heating 
and cooling terms sum up to zero in the equilibrium state: 

13 Actually this result only holds in the center of mass reference 
frame, i.e. with Vy Com ' subtracted. 
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-£c = 0. The work done by pressure forces 



balances dissipation by drag forces and losses from the 
backwards transport of angular momentum. 

5.1. Clumping and Dissipation 

In this subsection we will show that particle clumping 
reduces energy dissipation by drag forces, at least in a 
laminar state. Particles effectively "draft" off each other 
like birds flying in formation or bicycle riders in a pelo- 
ton. This drafting does not rely on overlapping turbulent 
wakes, but instead depends on slowing relative gas mo- 
tions by the collective inertia of particles. It is tempting 
to argue that the lowered dissipation rate explains the 
tendency of particles to clump. As usual, the story is 
more complicated, but the evolution of £drag turns out 
to be a useful diagnostic for the non-linear simulations 
of JY. 

First we demonstrate that dissipation is reduced by 
clumping. Consider the equilibrium drag dissipation of 
equation (|26|) . for simplicity in the tight coupling limit 
(t s <C 1), which we now express per unit surface area 
instead of volume as 



Adiss — ^diss-^p ~ 



4r s 



+ 



(30) 



Now imagine concentrating the particles into a volume 
smaller by a factor n > 1 via vertical setting or clumping. 
Compared to the uniform solids-to-gas ratio e the new 
value is ne in clumps and in voids. The new height- 
averaged dissipation rate is 



A* 

•'Miss 



4r B 



■Z p f2. 



(1 + ne) 2 

The fractional change in dissipation (for t s <C 1), 



A* 

y Miss 



diss 



l + e 
1+ne 



< 1 . 



(31) 



(32) 



shows that clumping decreases the net dissipation of well- 
coupled particles and that the effect becomes stronger 
with increasing e. 

Unfortunately there is no reason to expect in general 
that the dissipation rate decreases, especially since the 
system is not closed, but driven by pressure gradients. 
Examples of driven systems in which mechanical dissipa- 
tion increases with the spontaneous transition from lami- 
nar to turbulent flow include drag on a rigid body (e.g. an 
airplane wing) and Rayleigh convection with fixed tem- 
perature on the endplates (Jeremy Goodman, personal 
communication). Indeed the non- linear simulations of 
JY find that |£drag| could increase or decrease in the non- 
linear state. Obviously drag dissipation is affected not 
just by clumping (as in the toy laminar calculation here) 
but by the turbulent velocities that tend to increase dis- 
sipation. Nevertheless JY demonstrate that runs with 
the largest (and longest lived) overdensities show a de- 
crease in | f drag |) lending credence to the hypothesis that 
drafting can augment particles' ability to clump. 

6. DISCUSSION 



This paper begins our numerical exploration of the 
streaming instability, which uses aerodynamic particle- 
gas coupling to tap the radial pressure gradient in pro- 
toplanetary disks. Growing oscillations arise in an ideal- 
ized model for protoplanetary disks that assumes a local, 
unstratified, and non-self-gravitating shearing box with 
gas and uniformly-sized, non-colliding solids. Studying 
a relatively simple system isolates the surprisingly rich 
consequences of mutual drag coupling in disks. Also, 
the well-defined growth rates of seeded eigenvectors make 
the streaming instability an ideal test of numerical imple- 
mentations of particle-gas dynamics, as suggested in YG. 
We encourage those who study manifestations of particle- 
gas dynamics in disks to consider the linear streaming 
instability as a test problem if the feedback of solids on 
gas dynamics is relevant. 

This work is largely successful in reproducing the an- 
alytic growth rates of YG. The two-fluid simulations, 
which treat solids as a pressureless fluid, give excellent 
results with minimal computational effort. Particle- fluid 
simulations also converge to the analytic results, but 
higher spatial resolution is required. Treating the solids 
as particles has several advantages - it is more realis- 
tic, it can validate the often-used fluid approximation for 
solids, and it allows the development of non-linear den- 
sity enhancements without spurious shoc ks. Refinements 
of the particle-fluid algorithm used in fjohansc n et al.l 
(2006) are described, notably the use of higher order 
interpolation and assignment schemes to minimize er- 
rors in the drag force computation. These errors become 
more drastic as the stopping time decreases and errors 
grow relative to the diminishing difference between gas 
and particle velocities. Smaller stopping times also give 
shorter length scales, thereby imposing stricter Courant 
criteria. These restrictions actually dominate the obvi- 
ous concern that tighter coupling stiffens the equations 
of motion. Detailed modeling of the smallest particles 
in protoplanetary disks, especially in the inner regions 
with high gas densities, will require further algorithm 
development and increased computational power. In the 
meantime, studies of moderate coupling can establish the 
relevant physical phenomena and provide a baseline for 
extrapolation to more extreme parameters. 

Having developed a particle-mesh scheme that can be 
trusted to simulate coupled particle-gas dynamics with 
feedback, we proceed to explore the non-linear evolu- 
tion of streaming instabilities in a companion paper (Jo- 
hansen & Youdin 2007) with particular attention to the 
growth and saturation of particle overdensities. 
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APPENDIX 



WEIGHT FUNCTIONS AND INTERPOLATION/ASSIGNMENT ERRORS 



The choice of weight functions for interpolation, W\, and assignment, involves trading computational cost against 
performance. We will consider only 1-D weight functions which are combined multiplicatively to form multidimensional 
weights, W(x — xj) = W(x — xi)W(y — y m )W(z — z n ), for the cell centered at Xj = x^x + y m y + z m z. This assumes 
a rectilinear domain of influence, which is simpler (if less physical) than circular/spherical clouds. 

The interpolation function in the non- linear simulations uses quadratic splines (QS), 



(QS) 
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where Sxi = x — xi measures the distance from a cell center and Wj extends over three cells of width A. The 
interpolation errors are calculated by considering a periodic function (of arbitrary phase) sampled at the grid points. 
The QS interpolated values at arbitrary x in cell £ are 

(Afc) 2 i . Jl , r(Afc) 3 



cos(kx + <p) 
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The amplitude of a periodic signal is reduced by 1 — (Afc) 2 /8. For an arbitrary distribution, this smoothing can be 
simply corrected for in Fourier space. This sharpening is included, but only for the initial amplitudes of u, in some 
linear tests ( SI4.2|) . but was too costly for the non-linear runs in JY. There is also a noise, i.e. an error that depends on 
position relative to cell center §x#, of amplitude (/cA) 3 /(72v / 3). 

We also considered quadratic polynomial interpolation (QP), which performs a best fit through the three nearest 
grid points, resulting in a weight function 
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The QP interpolated values of a periodic signal read 
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The amplitude is preserved to second order, an improvement over quadratic spline interpolation. The noise, however, 
has an amplitude of (Afc) 3 /16, a factor of 15 larger than with quadratic spline. Also troubling is the discontinuity of 
(Afc) 3 /8 at the cell boundaries. 

The assignment function used in all the non-linear simulations of JY is the Triangular Shaped Cl oud (TSC) scheme, 
as opposed to the lower order NGP (Near est Grid Point, which was used in lJohansen et aL| [2006) or CIC (Cloud In 
Cell) schemes (|Hocknev &: Eastwoodlll98ll) . The TSC assignment weight function is identical to the quadratic spline 



interpolation weight function, W 



(TSC) 
A 



(QS) 



Assignment errors depend partly on sampling, i.e. the number of 



particles that make a non-zero contribution to a sum like equation (|19p . A higher order method like TSC samples more 
particles at a given grid point and gives a smoother distribution than lower order methods. The fractional amplitude 
reduction of a mode perfectly sampled by TSC is identical to the result for quadratic spline interpolation, 1 — (Afc) 2 /8. 
This is the same order, but larger than for CIC [1 — (Afc) 2 /12] or NGP [1 — (Afc) 2 /24] assignment, although especially 
the NGP scheme would require an enormous particle number to achieve good sampling of linear perturbations. In 
principle Fourier sharpening could be applied in the force assignment step, but it is less important to consider, since 
the errors are not magnified by a subsequent subtraction. 

NUMERICAL TEST OF PARTICLE ASSIGNMENT OVER SHEARING BOUNDARIES 

This Appendix studies the behavior of a linear, non-axisymmctric wave of gas and particles in order to test drag 
force assignment across the shear-periodic radial boundary. Solutions from the full simulation of the Pencil Code 
are compared to the following semi-analytic problem. We use local, linearized equations of continuity and motion to 
describe the evolution of gas and particle density p' s (x, y, t), p' p (x, y, t) relative to a constant density background state 
Pg,o, Pp,o, and gas and particle velocities u'(x,y,t), w'(x,y,t) relative to the Keplerian shear flow Vq = —(3/2)f2xij. 
We assume a shearing wave solution, q'(x,y,t) = q(t) exp\i( k x (t)x + k y y)], for each perturbation variable, with (see 
e.g. [Goldrcich fc Lvnd cn-Bell 1965; Brandenburg ct al. 200J) 



k x (t) = k x (0) + (3/2)ntk y . 
The wave amplitudes then evolve as coupled ordinary differential equations in time 14 , 



dpg/d< = -p g ,o [ik x (t)w x 



lkyUJy 
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Note that for this test problem no global pressure gradient, and thus no drift motions, are included. 
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Fig. 7. — The temporal evolution of a leading shear wave of gas and solid particles. The plot shows a comparison between the semi- 
analytic solution to the linearized equation system (gray lines) and the solution obtained with the full solver of the Pencil Code for the 
amplitude of the particle density p p (black dash-dotted line) and of the particle velocity components w x (black dotted line) and w y (black 
dashed line). There is excellent agreement between the numerical and the analytical solutions up until t ~ 5i? _1 where damping of the 
amplitude of the tightly wound wave by the Triangular Shaped Cloud scheme becomes significant. 



du x /db = 2Qu y - e (u x - w x )/t { ~ ik x (t)(%p & /p St0 , (B3) 

du v /dt = -Qu x /2 - e (u y - w v )/t { - ik y c^p g /p g ^ , (B4) 

dp p /dt = -p p ,o[ik x (t)w x + ikyWy] , (B5) 

dw x /dt = 2Qw y - (w x - u x )/t{ , (B6) 

dw y /dt = — Qw x /2 — (w y — Uy)/rt . (B7) 

where eo = P P .o/Pg,o- We solve this system of ordinary differential equations numerically for J? = p g ^Q = p P ,o = 
Tf = c s = k y = 1 using a third-order Runge-Kutta time integration method to follow the temporal evolution of a 
non-axisymmetric wave with the initial condition k x = — l,u x = p g = w x = w y = p p = 0,u y = 10~ 3 . The semi- 



analytic solution is then compared to the evolution obtained with the full solver of the Pencil Code using 64 2 grid 
points with 1 particle per grid point to cover a box of size L x = L z = 2ir. Fig. [7] shows the evolution of the absolute 
value of the particle amplitudes p p (dash-dotted line), w x (dotted line) and w y (dashed line) in comparison with 
the analytical solution (gray lines). There is an excellent agreement for t < 5.0. At later times, the wave becomes 
so tightly wound that damping of the wave amplitude by the TSC scheme becomes significant. Most importantly 
this non-axisymmetric test problem never shows any spurious features near the radial boundary (or anywhere else), 
validating our implementation of drag force assignment over the boundaries. 

COLD START: ALGORITHM FOR SEEDING DENSITY PERTURBATIONS 

Seeding low amplitude (we use 6 P = 10~ 6 ) density perturbations with particles is non-trivial. The desired density 
distribution cannot be seeded by random numbers for a reasonable number of particles, N p . The white Poisson noise 
has a constant Fourier amplitude of ~ 1/WN P at all scales, i.e., we would need a total number of particles N p S> 10 12 
to resolve 5 P — 1CP 6 ! 

Instead we borrow a tactic from cosmological simulations (e.g. iTrac fc Penll2006h to concentrate power in a desired 
mode. We first assign particles to a uniform grid with positions, Xi, labeled by a particle index i — 1,2, ...,N p . 
This grid is defined relative to the gas grid with an integer number of particles in each gas cell. We introduce linear 
perturbations to the density by applying periodic shifts to the particle positions. To approximate a density distribution 

p p (x) = (p p )[l + Acos(k ■ x)] (CI) 

with A <ti 1, the desired shift from the uniform grid is 

€< = -^Asm(k -x l ). (C2) 



14 



Youdin & Johansen 



The resulting density distribution, 



Pp( x ) = y2 s ( x - x i-€i), ( C3 ) 



has a Fourier transform 

Pp( k ) = VJi E ex P i ifc • [*< + &]} (C4) 



« [4,o + A(5 fe>feo + 4,- fco )/2] + p< 2 >(fc) + 0(A 3 ) , (C5) 

Vcell 

where V^ c u is the volume of a grid cell. The final step shows that we reproduce the desired plane wave to lowest 
order by performing an expansion about ^ <C Xi and using summation relations for periodic functions (we ignore 
the sub-gridscale aliases of ko). The standing wave solutions described in Sect. 12.31 are produced by summing two 
plane- wave displacements of ±k z . The quadratic error term 

M 4 2 

P ( p 2) (k) = ^ — (S k , 2ko +5 k ,- 2ko ) (C6) 
is small for A <C 1, but is eliminated by a further displacement 



Ci 2) = ^ 2 ^-'^„a-;.. (07) 



Thus equation (|C3|) has the desired Fourier properties even with only one particle per grid cell. However, the binned 
density distribution is actually what is relevant for influencing gas dynamics. When the TSC assignment scheme 
(described in Appendix B) is applied to the "cold start" positions, we cleanly get the desired density distribution 
assigned on the mesh, even for arbitrarily small shift amplitudes. 
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